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Abstract 

In  a  previous  paper,  we  introduced  a  coordinate-splitting  ( CS )  form  of 
the  equations  of  motion  for  multibody  systems  which  together  with  a  modified 
nonlinear  iteration  (CM),  is  particularly  effective  in  the  solution  of  certain 
nonlinear  highly  oscillatory  systems.  In  this  paper,  we  examine  the  convergence 
of  the  CS  and  CM  iterations  and  explain  the  improved  convergence  of  the  CM 
iteration.  An  example  is  given  from  flexible  body  simulation  which  illustrates 
the  convergence  results  and  the  class  of  problems  for  which  the  CM  iteration 
is  most  effective. 
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1  Introduction 


In  a  recent  development  for  the  solution  of  constrained  multibody  systems,  we  have 
proposed  a  coordinate-splitting  ( CS )  formulation  [24],  Numerical  experiments  have 
shown  that  the  CS  formulation  is  effective.  In  addition,  its  variant  CM ,  in  which 
the  second-order  derivative  of  the  C5-projection  operator  is  omitted  in  the  iteration 
matrix,  has  exhibited  surprisingly  good  results  for  the  solution  of  multibody  systems 
with  high-frequency  oscillations.  In  this  paper,  we  present  convergence  results  for 
the  CS  and  CM  method  for  the  discretized  system  of  equations  which  explain  the 
observed  behavior. 

First,  we  give  some  background  on  the  numerical  solution  of  constrained  multibody 
systems.  For  more  details,  we  refer  to  [2,  6,  8,  12,  17].  The  equations  of  motion  of  a 
constrained  multibody  system  can  be  written  as  [10] 

M(q)q-f(q,q,t)  +  G(q)T  A  =  0  (la) 

$(?)  =  0  (lb) 


where  q  =  [91,  q2i  9n]  316  the  generalized  coordinates ,  A  —  [Aj,  A2,  ...»  Am]  are  the 
Lagrange  multipliers ,  M(q)  €  JR"*"  is  the  mass-inertia  matrix,  /  €  SC  is  the  force 
applied  to  the  system,  q  =  %  is  the  velocity  and  q  =  0  is  the  acceleration  vector. 

The  constraints  g  —  [<71,  <72 >  •••■>  5m]  are  m  smooth  functions  of  <7,  whose  Jacobian 

G(q)  =  f|^l  €  Srxn  ,m<n  (2) 

oqj. 

is  assumed  of  full  row-rank.  We  also  assume  that  G(q)M(q)GT(q)  is  symmetric  and 
positive  definite  for  every  q  €  EC  to  obtain  a  consistent  physics  represented  by  (1). 

The  degrees  of  freedom  for  the  system  (1)  is  then  p,  where  p  —  n  —  m. 

It  is  useful  to  note  that,  for  a  conservative  system,  the  constraint  reaction  force 
GtX  in  (la)  is  obtained  by  the  derivative  of  reduced  potential  energy  [15] 

V*(«)  =  Ars(«)  (3) 

for  every  solution  of  A  and  q  from  (1).  Differentiating  Eq.  (lb)  yields  the  velocity  » 

constraints, 

G(q)v  =  0  (4) 

and  differentiating  Eq.  (4)  gives  the  acceleration  constraints, 

G(q)v  +  d^~v  s  G(q) v  -  y(v,  q)  =  0.  (5) 
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To  eliminate  the  constraints,  we  may  choose  m  coordinates  of  q  =  Xx  +  Yy  such 
that  (G(q)Y)~1  exists  in  a  neighborhood  of  q,  where  X  €  JR?Xp  and  Y  €  JR?xm, 
whose  columns  constitute  the  standard  basis  for  JR Since  G{q)  is  full-rank,  we  can 
use  the  splitting  of  q  —  Xx  +  Yy  to  obtain  a  linear  operator  defined  as  follows: 

Definition  1  [Coordinate-Splitting  Matrix]  Let  X  andY  be  the  matrices  whose  columns 
constitute  the  standard  basis  of  JR?Xn  such  that  ||(G,(g)y)“1||  is  bounded  in  a  neigh¬ 
borhood  Uo  of  go.  The  p  xn  coordinate-splitting  matrix  for  (1)  is  defined  by 

P(g)  =  XT-  Q(g)TYT  =  XT(I  -  G(q)T(G(q)Y)-TYT)  (6) 

where  Q(q )  =  {G(q)Y)-'G{q)X. 

Remark  1  From  the  construction  of  the  CS  matrix  P(q),  we  can  easily  see  that 
P(q)(jr(q)  =  0  for  all  q  €  JR?,  i.e.,  P(q)  is  orthogonal  to  range(GT).  Furthermore, 
the  row  vectors  of  P(q)  are  orthonormal,  i.e.,  P(q)TP(q)  =  Ip  where  Ip  is  the  identity 
matrix  in  JR?. 


Reducing  (1)  to  a  first-order  DAE,  appending  the  velocity  constraints  (4)  to  (1),  and 
applying  P(q)  to  its  differential  part,  we  obtain  an  index  1  DAE 


p(q)(q  -  v) 

=  0 

(7a) 

P(q)(M(q)v-f(v,q,t)) 

=  0 

(7b) 

G(q)v 

=  0 

(7c) 

=  0. 

(7d) 

Note  that  one  can  use  the  generalized  coordinate  partitioning  method  to  reduce  (1) 
to  the  p  differential  equations  [22] 

M(x,  h{x))x  =  f(x,  ^ x ,  x,  h(x),t)  (8) 

where  h(x)  is  the  implicit  function  of  y  defined  by  the  constraints.  The  equation  (8)  is 
equivalent  to  the  underlying  state-space  form  of  the  coordinate-splitting  formulation 
(7)  of  the  constrained  multibody  systems.  However,  the  function  h:JR?  — »  JR m  and 
its  time  derivative  cannot  be  evaluated  unless  the  constraint  equations  (lb)  and  their 
derivatives  have  been  satisfied,  e.g.,  the  computations  of  M  and  /  require  (x,h(x)) 
to  lie  on  the  constraint  manifold 

M  =  {q€JRn\g(q)  =  0} 
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and  %.*)  *  tangent  vector  to  M  at  (x,  h(x)).  These  conditions  necessitate  the 
solution  of  the  nonlinear  constraint  equations  and  their  derivatives.  The  numerical 
solution  of  (8)  can  be  inefficient  when  /  is  stiff ,  even  if  implicit  numerical  integration 
methods  are  applied,  because  of  the  computational  complexity  of  h(x),  /,  and  their 
derivatives. 

In  any  case,  it  is  advantageous  to  use  (7)  rather  than  (8)  for  the  numerical  solution 
of  (1)  because  there  is  less  computational  complexity  to  obtain  (7)  than  (8).  Applying 
numerical  integration  to  (7),  we  will  show  in  Sec.  2  that  the  CS  and  CM  iterations 
solve  the  nonlinear  system  (7)  efficiently  in  a  two-stage  iteration.  First,  we  carry 
out  the  iteration  of  a  2n  x  2n  system  of  the  Newton-Euler  equations  with  some 
additional  terms  corresponding  to  the  derivative  of  the  reduced  potential  (3).  Then, 
we  solve  for  the  increments  of  the  dependent  coordinate  y  holding  the  independent 
coordinate  fixed  on  the  range  space  of  the  projection  operator,  and  finally  update 
the  factorization  of  the  m  x  n  constraint  Jacobian  G.  For  the  CM  iteration,  the 
first  iteration  of  the  2n  x  2n  Newton-Euler  equations  (7a)  and  (7b)  uses  only  their 
unconstrained  form,  e.g.,  the  second-order  derivatives  are  omitted. 

In  Sec.  3,  we  show  the  convergence  of  the  CS  and  CM  iterations  under  a  mod¬ 
erate  assumption  on  the  smoothness  of  the  constraint  manifold  A4.  For  a  sufficient 
condition  of  the  convergence,  we  give  a  required  stepsize  of  the  numerical  integration 
methods.  The  sufficient  conditions  of  the  numerical  integration  for  the  CS  and  CM 
iterations  are  different,  given  the  same  starting  values  of  the  iterations.  For  problems 
with  a  small  potential  energy  force,  e.g.,  VVr(g)  =  is  small,  the  CM  iteration  is 
advantageous  over  the  CS  iteration.  As  shown  in  [24],  the  CM  iteration  has  not  only 
less  computational  complexity,  but  also  better  convergence  for  the  bushing  example, 
compared  to  the  C S  iteration. 

In  the  presence  of  highly  oscillatory  forces,  the  stepsize  h  of  the  numerical  inte¬ 
gration  may  be  restricted  to  h  <  y/\\^\\^  if  one  wiU  resolve  ^  the  high'frequenCy 
oscillations.  Here,  we  axe  not  interested  in  small  oscillations  with  high-frequency, 
hence  we  used  a  discretization  method  with  strong  damping  properties  and  the  step- 
size  h  is  not  restricted  by  the  high-frequency  oscillations  of  an  amplitude  smaller  than 
the  error  tolerance.  However,  the  convergence  of  the  Newton  iteration  for  (7)  is  sub¬ 
ject  to  a  good  approximate  Jacobian  at  the  predicted  solution.  Since  for  the  highly 
oscillatory  mechanical  systems,  the  Jacobian  evaluated  at  a  predicted  solution  may 
be  far  away  from  the  Jacobian  at  the  corresponding  numerical  solution,  the  Newton 
iteration  for  the  discretized  nonlinear  equations  has  convergence  difficulties.  For  some 
model  problems,  while  the  stepsize  in  the  CS  iteration  must  be  restricted  to  obtain 
Newton  convergence,  it  is  not  the  case  in  the  CM  iteration.  Using  a  simplified  New¬ 
ton  iteration,  we  analyze  the  rate  of  convergence  for  the  CS  and  CM  iterations  in 
Sec.  4  and  explain  the  different  convergence  behavior  for  the  highly  oscillatory  case. 
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In  modeling  a  deformable  body,  the  most  commonly  used  technique  is  the  finite 
element  method,  which  yields  linear  deformation  forces  in  the  body-fixed  local  coor¬ 
dinate  systems.  We  illustrate  that  the  theorem  in  Sec.  4  can  be  applied  to  this  class 
of  constrained  multibody  systems  and  predicts  the  results.  We  re-examine  the  2D 
bushing  problem  in  Sec.  5. 

2  Solving  the  nonlinear  system 

In  this  section,  we  examine  the  iterative  solution  of  the  C S  formulation  and  its  variant, 
the  CM  iteration.  Denoting  the  current  time  t  =  tn  and  (qn,  vn)  the  numerical 
solution,  applying  for  example  a  BDF  formula  to  (1)  yields 


P(Vn){phqn  -  vn)  =  0  (9a) 

P{qn){M(qn)phVn  —  f(vn,qn,tn))  =  0  (9b) 

G(9n)t>n  =  0  (9c) 

$(?n)  =  o  (9d) 


where  ph  is  the  discretization  operator.  We  will  investigate  Newton-type  methods  for 
the  solution  of  (9).  To  form  the  Jacobian,  we  will  need  to  find  the  derivative  of  a 
vector  function  in  the  form  of  P(q)r  with  respect  to  q.  In  [24],  it  was  shown  that  this 
can  be  written  as 

2j(P(  4)r)  =  with  »  =  (GU)Y)-rYTT.  (10) 

Using  the  product  rule,  we  obtain  the  derivative  of  P(q)r(q), 

t,™ <n> 

where  s  is  the  same  as  in  (10).  Thus,  for  a  given  CS  matrix  at  q ,  the  vector  function 
P(q)r(q)  is  differentiable  with  the  order  no  less  than  the  minimum  between  those  of 
r(q)  and  G(q). 

Using  (11),  and  denoting  ri  =  phqn  -  v„  and  r2  =  M(qn)phvn  -  f{vn,qn,in),  the 
Jacobian  of  (9)  is 
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where  Si  =  {GY)~TYTr \  and  sj  =  {GY)~tYtT2.  The  Newton  equations  of  the 
discretized  form  of  (9)  are 

P.  ([*£*  +  -  Ar„)  =  -P„r,  (13a) 

J>.  j[l|y)  +  =  -P„r2  (13b) 

^AS„  +  Gn&Vn  =  -G,Vn  (13C) 

GnAqn  =-9n  (13d) 


at  the  time  t  =  where  A qn  and  Av„  are  the  increments  of  qn  and  vn  by  the 
Newton  iterations.  For  notational  simplicity,  we  write  the  subscript  n  of  a  function 
representing  its  numerical  value  at  tn,  e.g.,  gn  —  g(qn ). 

A  modification  of  the  iteration  matrix  (12)  leads  to  the  CM  method  as  explained 
in  the  following.  Combining  (13a)  and  (13b),  we  obtain 

’o(,n)  V)](Ji(?"'“")[aS]  +  [S])=0  (14) 


where  the  2n  x  2n  matrix  J\  is 


J\{Qnivn)  — 


j£h3n  4.  rf(gn»i)  _/ 

dqn  '  dqn  1 

_  iZa  +  MdPhv"  _  ILl 

\n  dqn  dqn  dvn  dvn 


Replacing  by  P(q )^1,  i.e.,  fixing  the  CS  operator,  (15)  yields 


Ji(qn,vn)- 


_  M  dphVn  — 

dqn  dvn  dvn 


where  the  second-order  derivatives  and  in  (15)  are  nullified.  Replacing 
the  Jacobian  matrix  Ji{qn>  vn)  by  this  approximate  Jacobian  matrix  J\(qn,  t>n)  in  the 
iteration,  we  defined  the  CM  iteration  [24]. 

The  iterative  solutions  of  CS  or  CM  for  (9)  require  the  solution  of  the  linear 
system  (13),  which  can  be  obtained  from  one  additional  matrix  factorization,  i.e.,  the 
factorization  of  (15)  or  (16).  Since  J\  is  generally  invertible  under  the  assumption  of 
M(qn)  nonsingular,  one  solution  of  (14)  can  be  computed  by 

(XT  -  QlYT)(A,«  -  ft)  =  0  (17a) 

(XT  -  QlYT)(Av„  -  r2)  =  0  (17b) 


where 


Mqn,  <>n)  [  r\  ]  ”  [  -r\  ] 
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Note  that  (17a)  is  not  a  necessary  but  a  sufficient  condition  of  the  solution  of  (14). 
Then,  we  solve  (17a)  and  (13d)  for  A qn-  Using  A qn  to  compute  ^(?n)A qn,  we  solve 
(17b)  and  (13c)  for  Avn.  The  solution  (Ag„,  Avn)  requires  only  solving  two  linear 
systems  of  the  form 

[  G  ]  u  =  [  b  ]  (19) 

for  v  e  JR”,  a  €  .Rp,  and  b  €  JR”*.  Denoting  u  —  Xux  +  Yuy,  we  obtain  from  the  first 
p  equations  of  (19) 

ux  —  QTuy  =  a  (20) 

and  the  second  m  equations  yield 

uy  =  (GY)-\b-(GX)ux).  (21) 

Substituting  (21)  into  (20)  and  solving  for  ux  yields 

u,  =  (It  +  QTQ)-'(a  +  QT(GY)-'l')  (22) 

where  Ip  is  the  p  x  p  identity  matrix.  According  to  (21)  and  (22),  the  solution  of  (19) 


depends  only  on  its  independent  part,  e.g.,  ux. 

Combining  (17a)  with  (13d),  we  obtain 

A*.  =  -(/,  +  +  Ql^n  -  (GnY)~lgn))  (23a) 

Ay„  =  -(G„y)"1(jn  +  (G„X)Ai„)  (23b) 

and  from  (17b)  and  (13c),  we  have 

Au>n  =  -(/P  +  Q^n)-1(-^rr2  +  Ql(Frr2-(Gny)-1»?„))  (24a) 

Azn  =  -(GnY)-\r,n  +  (GnX)Awn)  (24b) 


where  rjn  =  Aqn,  and  vn  =  Xwn  +  Yzn.  The  numerical  solutions  (23)  and  (24) 
illustrate  that  the  dependent  variables  y„  and  zn  are  determined  geometrically,  e.g., 
use  only  the  algebraic  constraints.  Therefore,  applying  numerical  integration  to  qn 
and  vn,  the  local  errors  can  be  bounded  by  the  difference  in  xn  and  wn  using  CS  or 
CM  iterations. 

i  It  is  useful  to  examine  the  difference  between  the  CS  and  CM  iterations.  The 

components  of  J\  that  are  dismissed  by  the  CM  iteration  represent  the  derivative  of 
»  the  CS  matrix  P(q),  that  is,  the  tensor  of  the  second-order  derivative  of  the  constraint 

equations.  The  influence  of  this  tensor  on  the  increment  Aq  may  be  expressed  by 

l(P(,)r)A,  =  r  (25) 
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where  Ay  =  (GY)~TYT Aq.  The  exchange  of  Aq  and  r  in  (25)  is  permitted  by  the 
smoothness  of  the  constraint  manifold  defined  by  M  =  {9  €  lR^\g(q)  =  0},  on  which 
is  a  bilinear  form.  The  term  (25)  measures  the  rate  of  change  of 

the  normal  vector  C?Ay  (to  the  constraint  manifold  at  q)  along  the  solution  curve 
on  the  independent  generalized  coordinate  space.  The  difference  between  the  CS 
and  CM  iterations  can  be  expressed  in  terms  of  (25)  with  the  corresponding  residual 
vectors  r\  and  r2  of  (13a)  and  (13b),  respectively. 


3  Convergence  of  CS  and  CM  Iterations 

The  convergence  results  for  the  CS  and  CM  iterations  can  be  carried  out  on  a  smooth 
constraint  manifold  «M.  We  assume  that  for  any  qo  €  At ,  there  exist  X  €  lRpxn  and 
Y  €  ]RmXn  such  that 

\\(G(q)Yr\\<C1  (26) 

||G(9l)T  -  G(?2)T||  <  C2\\qi  -  fe||  (27) 

for  some  C\  and  C2,  where  q,  91,  and  92  are  in  a  neighborhood  U(qo)  of  90- 

Remark  2  The  matrices  X  and  Y  for  a  given  q0  may  be  selected  according  to  dif¬ 
ferent  strategies.  For  instance,  applying  Gaussian  elimination  with  row  pivoting  to 
GT(qo),  one  obtains  Y  by  the  permutation  indices  of  the  factorized  matrix.  In  this 
case,  C\  of  (26)  is  the  same  order  of  magnitude  as  ||(GfGr)“1||.  On  the  other  hand, 
C2  of  (27)  is  the  Lipschitz  constant  of  G,  which  is  independent  of  the  choice  of  X 
and  Y. 


With  the  conditions  (26)  and  (27),  it  is  easy  to  obtain  an  upper  bound  for  the 
difference  term  (25)  between  the  CS  and  the  CM  iterations. 

Lemma  1  Suppose  conditions  (26)  and  (27)  hold.  Then 

l4[F(?)r(«)]-i>(«)^||<ft,C.CJ||KTr(9)ll  (28) 

in  D(q0,Ro)  Q  U(q0),  where  D(q0,  Ro)  is  the  disc  in  JR?  with  center  q0  and  radius 
Rq. 


Proof.  The  inequality  is  a  direct  consequence  of  (10)  and  (11).  Subtracting  (10) 
from  (11)  and  taking  the  norm  of  the  remainder  yields 
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Since  the  row  vectors  of  P(q)  are  p  orthonormal  vectors  in  IR*1,  applying  the  Cauchy 
inequality  gives 

|f(,)«l  <  ll«ll  <  RoCMGY)-TYTr(<,)\\ 
for  all  q  €  D(qo,Ro)  C  U(qo).  Condition  (26)  implies  the  result  in  (28).  □ 

For  simplicity  we  now  consider,  instead  of  the  second-order  constrained  equations 
of  motion  (1),  a  first-order  system  of 

«-/(«,<) +  GrA  =  0  (29a) 

,(,)  =  0  (29b) 


1 


% 


since  the  convergence  results  of  CS  and  CM  methods  for  (29)  can  be  trivially  ex¬ 
tended  to  (1).  Applying  coordinate-splitting  to  (29)  at  the  time  t  =  we  obtain  the 
nonlinear  equations 

W)  =  [  ]  (30> 

where  the  residual  function  is 

r(q,t)  =  Ph(q)  -  (31) 

using  the  linear  discretization  operator  ph  with  the  time  step  h.  Denoting  the  nu¬ 
merical  solution  of  the  nonlinear  equations  (30)  by  q *,  and  the  numerical  solution  of 
(31)  by  q*  at  the  time  t  =  t*  we  apply  the  theorem  of  Newton-Kantorovich  to  the  CS 
iteration,  yielding  the  convergence  of  CS  [5]. 


Theorem  1  (Convergence  of  CS  iteration)  Suppose  conditions  (26)  and  (27)  hold 
at  the  solution  q*  of  (SO)  for  a  selected  Y.  Applying  multistep  numerical  integration, 
if  the  numerical  discretization  operator  p^  satisfies 


<Co<l 


(32) 


with  s  =  (GY)~TYTr,  for  all  q  €  U(q*),  where  the  numerical  solution  q*  of  (SI)  is 
sufficiently  close  to  q*t  e.g.,  q*  €  U(qm)  H  U(qm),  then  for  all  starting  values  qo  such 
that  ||<7*  —  <7o||  <  Ro  for  some  Ro>0,  the  sequence  {g*}  generated  by 


qk+ i  =  qk-  J(qk)  F(qk)  =  qk- 


G(ft) 


1-1 


F(qk)  (33) 


where  r,  =  ^  and  a*  =  ( G(qk)Y)~TYTr(qk ),  converges  to  q \ 
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Proof.  Under  the  assumptions,  the  Jacobian  matrix  J (<?*)  of  the  equations  (30) 
is  nonsingular,  G(q)  is  Lipschitz  continuous,  and  ||(^(r  +  G7 5))-1||  <  Cq  <  1  in 
U(q*)  OU(q*).  Thus  J(q)  is  Lipschitz  continuous  and  J(q*)~1  is  bounded.  Therefore, 
{qk}  defined  by  (33)  will  converge  if  the  initial  guess  q0  is  sufficiently  close  to  q*.  This 
is  a  standard  proof  of  convergence  of  the  Newton  method  for  (30).  For  details  we 
refer  to  [5],  pp.  90*91.  □. 

Remark  3  The  assumption  q *  €  U(q*)  fl  U(q*)  and  the  existence  of  J-1  depend  on 
the  time  steps  h,  of  (SI).  Since  the  time  step  is  not  the  focus  here,  we  will  assume 
an  appropriate  h  for  pk  in  all  the  discussions.  In  general,  we  may  obtain 

||,-  -  ,-||  <  min{C„||r(,-)||,  C,||j(f)||}  <  0(A>) 

for  some  j  consistent  with  the  order  of  discretization  operator  pk,  since 

Ilf  -«(f)ll  <  0(h‘)  and  ||f  -  ,(C)||  <  0(h>) 

where  q(t*)  is  the  analytical  solution  of  (29).  The  Rq  in  Theorem  1  may  be  taken  to 

be  i 

ii°“2||.7->||||^||' 

For  convergence,  one  of  the  sufficient  conditions  for  the  above  theorem  requires 
that  the  numerical  integration  satisfies  (32).  This  implies  that  the  stepsize  h  has  an 
upper  bound.  For  linear  multistep  integration,  e.g.,  f ,  the  stepsize  must 

satisfy 

(34) 

where  a  is  the  leading  coefficient  of  the  numerical  integration  formula.  Convergence  of 
the  CM  iteration  can  also  be  assured  for  a  sufficiently  accurate  initial  guess.  Carrying 
out  the  CS  iteration  (33),  the  increment  A qk  €  SC  satisfies 

PM  =  -?(» M?*)  (35a) 

G(qk)Aqk  =  -g{qk)  (35b) 

at  the  kth  iteration.  The  corresponding  CM  solution  at  the  k th  iteration  yields 

P(<lk)r9(qk)Aqk  =  -P(qk)r(qk)  (36a) 

G(qk)Aqk  =  -g(qk).  (36b) 

A  bound  on  the  difference  between  Aqk  and  Aqk  can  be  computed  and  convergence 
of  the  CM  iteration  can  be  shown  as  follows. 
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(37) 


ft 


Lemma  2  Suppose  conditions  (26),  (27)  and 

hold  in  a  neighborhood  of  q*,  and 

^•)=[i’W)] 

is  invertible .  Let  {<?*}  generated  by  the  CS  iterations  converge  to  q*,  and  Aqk  be 
the  solution  of  the  CM  iteration  at  qk  for  k  =  0,1,2,....  Then  {6qk  =  Aqk  —  A g*} 
satisfies 

IM  <  ftC1C1||J-1(?-)lll|},rr(«,)lll|A?*ll  =  0{h’)  (38) 

using  a  jth-order  numerical  integration  method,  for  some  k  >  K\  >  0,  for  K\  suffi¬ 
ciently  large. 


Proof.  Subtracting  (35)  from  (36)  yields 


where 


This  implies 


ll%ll  <  l|J-1(ft)lll|ffe)<iG(’*)T(^)':rrTrt|ll|Aftll- 

Using  (28),  for  any  k  >  Ko,  the  first  term  of  the  right-hand  side  of  (38)  is  obtained 
by  (28).  To  show  that  this  term  is  actually  0(hJ),  we  note  that  ||A^||  =  0(h*)  for  k 
sufficiently  large,  since  \\Fk\\  — ►  0  when  k  —*  oo.  □ 


Theorem  2  (Convergence  of  CM  iteration)  Under  the  conditions  in  Lemma  2, 
choosing  qo  =  qo,  the  sequence  {qk}  generated  by  the  CM  iterations 

qk+i^qk-Jiqky'Fiqk)  (39) 


converges  to  q*  and 


\\r-q(n\\<o(v). 


(40) 
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Proof.  Since  J  is  nonsingular  and  its  components  are  smooth  functions,  using  (21) 
and  (22),  we  can  write 


for  the  CM  iteration,  providing  J  is  invertible.  By  conditions  (26)  and  (27),  we  have 


and 


for  some  constants  Cz  and 
gence  of  the  CM  iteration. 


h  QT{GY);'  1 1|  <  c4 
o  (<7y)-\  J11-  4 

Thus,  the  contractive  condition  (37)  implies  conver- 


□ 


Note  that  (37)  for  the  CM  method  is  analogous  to  (32)  for  the  CS  method.  Instead  of 
(34)  for  multistep  integration  methods  using  the  CS  iteration,  the  stepsize  condition 
for  the  CM  iteration  is  given  by 


h 

a 


(41) 


in  accordance  with  (37). 


Remark  4  The  sufficient  condition  (41)  is  not  a  necessary  conditioner  the  conver¬ 
gence  of  the  CM  iteration.  If  the  CM  iteration  is  carried  out  by  the  solution  of  (18) 
followed  by  (19),  then  the  stepsize  condition  (41)  can  be  rewritten  as 


It  is  easy  to  see  that  {qk}  of  the  CS  iterations  and  {qk}  of  the  CM  iterations  are 
the  same  if  g  is  linear  in  q.  In  general,  the  rate  of  convergence  of  the  CM  iteration  is 
superlinear ,  using  the  Dennis-More  Characterization  Theorem  [7].  Moreover,  if  the 
constraints  are  actually  invariant  to  the  differential  equations,  i.e.,  lim/,_o  r{q  )  =  0> 
then  we  have  J (q*)  —  J (?*)• 

It  is  noteworthy  that  the  CS  iteration  can  also  be  implemented  in  the  stabilized 
DAE  formulations,  which  are  based  on  the  application  of  the  method  of  Lagrange 
multipliers  [24].  For  example,  we  have  considered  the  stabilized  index-2  form  [8].  The 
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algebraic  variables  are  obtained  from  the  solution  of  si  and  in  the  computation  of 
the  Jacobian  (12).  Analogous  to  the  CM  modification,  we  can  eliminate  the  second- 
order  derivative  of  the  reaction  forces,  e.g.,  CF A  in  (lb),  yielding  an  approximate 
Jacobian  similar  to  that  of  the  CM  iteration.  However,  additional  modifications  to  the 
convergence  test  of  the  Newton  iteration  for  the  Lagrange  multiplier  based  formulation 
are  required.  For  the  Newton  convergence,  a  similar  result  to  Theorem  4  holds  for 
the  Lagrangian  formulations  in  certain  problems,  but  the  modified  convergence  test 
is  no  longer  reliable. 

4  Rate  of  convergence  for  highly  oscillatory  multi¬ 
body  systems 

One  challenging  class  of  problems  in  multibody  dynamic  systems  (MBS)  is  the  solu¬ 
tion  of  systems  with  high-frequency  vibration.  High  frequency  oscillatory  forces  often 
appear  in  the  modeling  of  vehicle  suspension  systems,  modal  analysis  in  structural 
dynamics,  or  modeling  oscillations  in  computer-aided  engineering  etc.  Using  the  CS 
formulation,  we  may  write  the  oscillatory  equations  of  motion  as 


PM(i-v)  =  o 

(42a) 

P(q)(M(q)v  +  -  /(«,?,<))  =  0 

(42b) 

o 

II 

(42c) 

o 

II 

CD 

(42d) 

where  ^  may  be,  for  example,  the  coefficients  of  stiff  springs;  i.e.,  0  <  £  <  1.  In 
practice,  q{q)  is  usually  oblique  towards  KerP(q),  i.e.,  the  oscillatory  force(s)  act  on 
both  the  independent  and  the  dependent  coordinates.  For  the  numerical  solution  of 
(42),  experiments  have  shown  that  the  CM  iteration  performed  superior  to  the  CS 
iteration  [24]  for  those  types  of  problems.  For  the  purpose  of  obtaining  a  smooth 
solution  with  large  stepsizes  [18],  we  can  explain  the  reason  why  the  CM  iteration  is 
so  effective. 

In  the  modeling  of  deformable  multibody  systems,  the  nonlinear  oscillatory  forces 
in  (42b)  are  usually  derived  from  the  theory  of  linear  elasticity,  i.e.,  for  some  functions 
q  such  that  the  oscillatory  forces  may  be  written  as  ^q.  We  can  use  these  functions 
q  to  write  the  nonlinear  force,  e.g. 

?»(«>  = 
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and  then  append 


q-Tj(q)  =  0 

to  the  constraint  equations.  The  oscillatory  forces  will  then  become  linear  with 
respect  to  the  variables  q.  In  fact,  if  the  oscillatory  forces  were  produced  by  a  finite 
element  approximation  of  the  deformation  of  bodies,  components  of  q  are  associated 
with  some  body-fixed  local  coordinates  via  the  orientation  transformation  matrix, 
whose  entries  often  are  slowly  varying  in  time. 


Deformation  forces  are  the  most  common  potential  forces  that  can  produce  small 
amplitude  high-frequency  oscillations,  and  they  are  usually  linear  with  respect  to  the 
local  coordinates  [4,  25].  For  these  reasons,  we  will  consider  the  class  of  oscillatory 
forces  in  the  form 

i(«) -«(*)(« -MO)  («) 

where  components  of  B  and  bo  are  slowly  varying.  In  particular,  B  and  bo  may  be 
functions  of  some  constraint-driven  generalized  coordinates.  For  example,  B(6)  in 
the  2D  bushing  problem  in  [24]  has  the  form 


B(0)  = 


COS0 

—  sin# 

0 


sin  6 
cos  6 
0 


kx 

0 

o  1 

r 

0 

** 

0 

0 

0 

k9 

. 

cos  6 
—  sin# 
0 


sin0  0 
cos  6  0 
0  1 


(44) 


where  6  is  small,  and  ib*,  Jbv  and  k9  are  positive  constants. 

Using  a  linear  oscillatory  force,  the  Lagrange  equations  of  motion  of  the  MBS  can 
be  written  as 

M(q)v  +  ±B(q  -b0)  +  GfT\-  /(v,  9,  t)  =  0  (45) 

where  ~  >>  lla^jll*  From  the  assumption  (27)  of  the  constraint  manifold,  we  can 
also  see  that 


-  >  max 
e  IMI.IMI=i 


.dG{q)ui 
1  dq 


«2|| 


(46) 


for  all  q. 


In  the  context  of  the  CS  formulation,  the  problem  of  convergence  of  the  Newton 
iteration  can  be  explained  by  analyzing  the  reduced  potential  function.  The  reduced 
potential  of  (42b)  is 

Vr(q)  =  g(q)T(GY)-TYTr  (47) 

where  r  =  /  —  Mq  -  \B(q  -  bo).  The  reduced  potential  force  generated  by  (47)  is 


VVr(q)  =  =  GT{GY)~TYTr. 

dq 


(48) 
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At  each  iteration,  the  reduced  potential  force  acts  along  the  normal  direction  of 
the  constrained  manifold  enforcing  the  constraint  equations.  The  gradient  of  the 
correction  term  yields 

V2Vr(q)  =  ( J  -  Cfr(GY)-TYT)^^  (49) 


where  s  =  (GY)~TYTr.  Applying  YT  to  (49),  gives 

YTV2Vr{q)  =  Yt(I  -  GT{GY)-TYT)dGT^q--  =  0 


and  applying  XT  to  (49)  yields 

XTV2Vr(q)  =  P(q)d-}—- 

When  high-frequency  oscillations  appear  in  the  system,  e.g.,  e  — *  0,  the  reduced 
potential  force  also  becomes  oscillatory  if  YTr  is  nonzero.  This  is  the  genera]  case 
when  the  solution  is  not  at  an  equilibrium  position.  Nevertheless,  convergence  of  the 
CS  iteration  can  be  achieved  by  using  a  small  enough  stepsize,  e.g.,  h  «  y/e. 

Theorem  3  Let  (q,v)  be  the  solution  of  the  nonlinear  system  (42),  which  results 
from  numerical  integration  using  pn  with  a  stepsize  h.  Suppose  the  starting  value 
(qo,v0)  satisfies  ||q0||  =  0(h2)  and  ||v0||  =  0{h),  and  J(q0)  is  nonsingular.  Then  the 
CS  iteration  converges  if  h?  <  ce  for  some  moderate  c. 

Proof.  For  the  convergence  of  the  CS  iteration,  we  need  to  show  that  (32)  is  valid, 
where  r(q)  is  defined  in  (45).  For  (32),  we  have 

udr,  s  ,  dGFs,  Nl)  „or||Af||  1„ 

l|ai(«,)  +  _5r(,o)ll  =  ll_A3 — -J  +  °W 

where  a  >  0  is  the  leading  coefficient  of  pn ,  and  ||M  ||  is  not  zero.  Consequently,  for 
e<l,  (32)  is  valid  provided  that  h  «  y/e.  □ 

From  the  above  theorem,  we  can  obtain  the  same  convergence  result  for  the  CM 
iteration  using  Lemma  2,  provided  J(q0)  is  invertible.  In  many  applications,  following 
the  oscillations  is  not  of  interest.  Instead,  one  wants  to  use  a  large  time  step  to  damp 
out  the  oscillations  of  small  amplitude  but  high  frequency.  For  this  reason,  we  now 
consider  only  the  multistep  numerical  integration  methods  that  axe  strictly  stable 
at  infinity  and  A-stable,  such  as  the  lower  order  (i.e.,  <  2)  BDF  methods  [12].  The 
convergence  of  Zr-stable  implicit  Runge-Kutta  methods  to  the  smooth  solution  of 
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highly  oscillatory  ODE  of  multibody  mechanical  systems  can  be  found  in  [18].  Here 
we  focus  on  the  convergence  of  the  C M  iteration  for  constrained  multibody  systems 
with  oscillatory  forces  when  applying  the  above-mentioned  linear  multistep  methods. 

Numerical  solutions  on  the  slow  manifold  can  be  evaluated  using  the  equilibrium 
of  (42b),  i.e.,  the  slow  solution  [1, 14]  satisfies 

V(q)  ~  -  Wr(g)  -  M(q)v)  =  0, 

and  the  smooth  solution  is  its  asymptotic  expansion  to  some  order  of  e  around  the 
manifold  {q  \  f](q)  -  0}.  In  the  linear  form,  the  smooth  solution  of  (42)  is  not  far 
from  B(q  -bo)  =  0  since  *  >  ||  For  the  strongly  damped  numerical  solution 

qn,  B(qn  -  bo)  ->  0(e)  as  tn  ->  oo.  During  the  iterative  solution  onto  the  slow 
manifold ,  the  constraints  may  not  be  satisfied,  which  causes  a  large  reaction  force 
in  the  form  of  (48).  This  may  cause  oscillations  in  the  CS  iteration,  while  the  CM 
iteration  annihilates  these  nonlinear  oscillations  generated  by  the  reduced  potential. 
This  yields  a  superior  performance  of  the  CM  iteration  as  compared  to  the  CS 
iteration  for  computing  the  smooth  solution  of  (42).  The  result  is  explained  in  the 
following. 


Lemma  3  Let  (qm,  u*)  be  the  smooth  solution  of  (42),  r}(q)  linear,  and  h  the  stepsize 
of  the  multistep  integration  method.  Suppose  the  starting  values  (q0,  u0)  for  (qm,vm) 
on  the  smooth  solution  of  (42),  i.e.,  ||9*||  =  0(c)  and  r(q\v’)  =  0(h),  satisfy  (26), 
(27)  and 

M(q0)ph(v0)  -  /( t>o>  9o)  =  0(h)  (50) 

where,  ph,  is  the  corresponding  discretization  operator.  Applying  the  CS  and  CM 
iterations  to  (48),  the  approximate  Jacobian  matrix  for  the  CS  iteration  satisfies 

l|J(flb.»>o)  -  =  |0(A)  +  O(A)  (51) 

where  S  =  \\Bq0  -  Bq*\\,  and  J(q,v)  is  the  Jacobian  of  (42).  For  the  CM  iteration, 

we  have  ,  . 

P(q°,vls)-J(q',v‘)\\=0(h)  (52) 

where  J  is  the  approximate  Jacobian  in  the  CM  iteration. 

Proof.  The  difference  between  the  Jacobian  at  (g0,t>0)  and  (<?*,u*)  can  be  written  as 

||Jo  -i’ll  <  |P(»)|(»)  -  W)fj(«”)  II  +  II^(9oM«o)  -  ^(9’M?-)II  +  °(A) 
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since  the  initial  values  satisfy  (50).  Under  the  conditions  (26)  and  (27),  we  may  choose 
common  X  and  Y  for  P(qo)  and  P(q*)  such  that  the  first  term  on  the  right-hand  side 
of  the  above  inequality  can  be  rewritten  as 

m?x^(9o)  -  ^(«-))ii = ow 

for  some  q  €  [?o> ?*],  since  +  0(h)  allowing  the  cancellation  of  ^B.  The 

second  term  yields 

II^(9oM«o)  -  ^(<I*)r(®-)ll  <  ||r(go)  -  r(,*)||0(A)  =  -  £5-||0(M 

according  to  Lemma  1.  Thus,  (51)  is  proved.  Recalling  J(qo,Vo)  from  (16),  we  have 

iUo-j-ii<imin^ii=ow, 

using  again  Lemma  1.  □ 


Theorem  4  For  the  initial  values  (qa,  t?o),  suppose  the  conditions  in  Lemma  3  hold. 
Suppose  that  the  CS  and  the  CM  iterations  are  carried  out  by  applying  a  simplified 
Newton  method ,  where  the  iteration  matrix  is  computed  at  the  starting  values  (qo,v 0). 
If  both  iterations  converge,  then  the  rate  of  convergence  of  the  CS  iteration  a(cs) 
compared  to  that  of  the  CM  iteration  <riCAf)  is  given  by 

alcs)  =  lo(h)  +  0(h) 

where  6  =  \\B(q0  — qm)\\,  and 

a(CM)  =  0(hy 


Proof.  We  consider  the  rate  of  convergence  that  is  defined  by 


a  —  lim  sup 
*-►00 


lk*+i-g*l 
Ik*  -  ?i  ' 


(53) 


Since  we  apply  the 
be  written  as 


where 


simplified  Newton  iteration,  the  solution  of  the  CS  iteration  can 
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Similarly,  the  CM  iteration  can  be  written  as  the  fixed-point  iteration  of  the  function 

H(q,v)  =  ]  “  JolF(q,v). 

Applying  the  Contractive  Mapping  Theorem,  see  [5]  pp.  93-94,  we  obtain  the  rates 
of  convergence  of  the  CS  and  CM  iterations: 

<rC5  =  P  -  tfjtfyy II  =  II  Jo 1 Vo  -  r)  II 

and 

aCM  =  ||J  -  J„-V(4-,V)||  =  II  Wo  -  -nil 

where  J(q*,vm)  =  J*  is  the  Jacobian  at  the  solution  of  the  discretized  system,  and 
the  superscripts  denote  the  respective  iterations.  For  Jq1,  we  have 

\  &  +  0(h)  -I  1 
J°  ~  [  \B  0(h)  l M  +  0(h)  * 

When  €  -+  0,  the  dominant  components  of  Jq1  are  of  0(1).  From  Lemma  3,  the  rates 
are  . 

acs  _  -^0(h)  +  0(h) 

and 

aC5  =  0(h) , 

since  Jq1  has  no  component  of  0(\).  □ 

In  the  stiff  bushing  example  of  [24],  we  have  seen  that  the  CM  iteration  was  far 
more  efficient  than  the  CS  iteration.  The  results  match  the  prediction  of  Theorem 
4,  e.g.,  as  e  — ►  0,  the  CS  iteration  became  very  ineffective  due  to  failures  in  the 
convergence  test  for  the  Newton  iterative  solutions. 

5  Example 

An  important  class  of  applications  of  the  CS  and  CM  iterations  is  flexible  multibody 
dynamics,  which  leads  to  the  coupled  large  displacement-small  deformation  equations 
of  motion  [25,  16].  As  shown  schematically  in  Fig.  1,  the  deformation  force  between 
the  i th  and  j th  components  is  a  function  of  the  relative  displacement  of  the  reference 
frames  X[-Yl-Z\  and  X'^Yj-Z'y  Typically,  the  relative  displacement  is  measured  by 

dij  =  rj  +  Ajs'j  -  r,  -  A, -s'  (54) 
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where  s[  and  Sj  are  constant  vectors  to  the  origins  of  the  force  reference  frames  in 
their  respective  local  coordinate  systems,  i.e.,  Xi-Yi-Zi  and  Xj-Yj-Zj ,  where  r,-,  rj 
are  the  corresponding  origins  in  a  global  coordinate  system  and  A{  and  Aj  are  the 
transformation  matrices  from  the  global  to  the  local  coordinate  system  [10].  The 
relative  angles,  6,j  =  Oij,  <f>ij]T,  are  calculated  as 

An  =  (AiBi)T  AjBj 

ipij  —  -A,j(2,3) 

^•  =  ^•(1,3) 

^•-arctan(^(2  2)) 

where  is  the  component  of  the  k th  row  and  /th  column  of  A(j.  The  matrix 

Aij  is  the  relative  orientation  matrix  of  two  force  reference  frames,  i.e.,  £,  and  Bj 
are  constant.  The  relative  velocity  is  the  time  derivative  of  the  relative  displacement 
dij  =  -fidij  and  relative  angular  velocity  is  =  u>j  —  o>,-,  where  u>i  and  Uj  axe  angular 
velocities  of  body  i  and  j,  respectively. 


Figure  1:  Deformation  Force  of  a  Flexible  Body 


Using  the  above  defined  notation,  the  force  acting  between  the  i th  and  yth  components 
due  to  the  deformation  can  be  written  as 


fii  =  AiB^K'iAiBifdu  +  C^AiBifiij) 
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where  Kf  is  a  3  x  3  structural  stiffness  matrix  and  C]  is  the  3  x  3  damping  coefficient 
matrix.  Similarly,  the  torque  acting  between  the  components  is 

Tn  =  AiBi(KTQij  +  C^AiBifua) 

where  Kr  and  Cr  are  analogous  to  Ks  and  Cf.  Note  that  the  force  and  torque  in 
this  form  are  linear  functions  of  the  relative  displacement  (d,j,  ©,j)  and  the  relative 
velocity  (d,j,fl,j). 

To  illustrate  the  deformation  forces  acting  on  bodies,  we  consider,  in  Cartesian 
coordinates  (x,  y,  0),  a  simplified  2D  bushing  force  at  the  body-fixed  coordinate  frame, 
whose  origin  is  (-5, 0)  and  the  axes  are  parallel  to  the  body-fixed  centroid  frame.  The 
other  reference  coordinate  frame  is  fixed  to  the  global  position  (■5, 0),  whose  axes  are 
parallel  to  the  global  coordinate  frame.  The  bushing  force  has  stiffness  matrix 

,  r  **  o  0 

Kf  =  -  0  &  0 

«  L  0  lb* 

where  kx,  kv,  and  ke  are  0(1),  and  damping  matrix 

c*  0  0  1 

Cf=  0  c*  0 

0  0  ce  m 

where  c*,  c*,  and  c9  are  0(1).  For  the  2D  bushing  example,  (44)  becomes 

.  r  1  0  0  "I  [  cos2  6  cos  8  sin  6  0  " 

B(0)  =  -  0  10  Ifc*  —  cos 8 sin 8  sin20  0  • 

e  [  0  0  lb*  J  [  0  0  0  J 

The  kinematic  constraints  on  the  body  are 

i(x2  +  k1  —  r)  =  0 

and 

8-8 o  =  0 

where  V  is  a  constant.  Let  the  gravity  and  the  mass-inertia  be  unity,  and  fc*  =  ky  =  1, 
then  the  equilibrium  of  the  bushing  force  satisfies 

x  =  ^(1  +  COS0O) 
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along  with  the  constraint  equations.  For  the  equilibrium  of  this  system  to  lie  on  the 
constraint  manifold,  we  should  have 

/*  =  e(e  — sin0o)  + ^(1 +cos0o)- 

For  simplicity,  we  consider  the  case  0O  =  0  and  the  corresponding  le  such  that  the 
equilibrium  (l,-e,0)  of  the  bushing  force  lies  on  the  constraint  manifold.  Near  the 
equilibrium,  we  can  choose  the  C S  matrix  as 


since  x  «  1.  Applying  a  numerical  integration  formula  (such  as  BDF)  to  the  CS 


formulation  of  the  system  yields 

phy  -w  -  -{phX  -  v)  =  0  (55a) 

X 

PhW-fy  +  l--(PhV-fx)  =  0  (55b) 

X 

xv  +  yw  —  0  (55c) 

u>  =  0  (55d) 

i(x2  +  y2-/')  =  0  (55e) 

6  =  0  (55f) 

where  the  velocity  is  [t>,u>,u>],  and  the  bushing  force  is 

'  /*  1  r  \  -x  +  \cnsO  ' 

fy  =  B(6)  -y  +  isin0  . 

f  -V  J 


Using  the  starting  values  x0  =  0.9,  y0  =  -yfi  ~  *?,  the  velocity  [-4.843  x 

10"3,  -1.0  x  10-2,0],  we  obtain  the  results  of  the  CS  and  CM  iterations  from  DASSL 
[19]  with  various  initial  stepsize  ho  and  the  ending  time  at  ho.  The  solution  tolerance 
TOL  is  set  to  10-10,  and  the  stiffness  e  =  10"5  and  the  damping  factors  (?  =  (?  =  10. 
As  shown  in  Table  1,  using  only  the  first-order  Euler  method,  the  CS  and  CM 
iterations  are  nearly  identical  for  the  initial  stepsize  h0  ranging  from  10-3  to  0.05. 
This  is  because  6  is  small  due  to  the  fact  that  the  predictor  is  very  accurate.  The 
approximate  rates  of  convergence  of  the  two  iterations  when  read  in  DASSL  agree 
up  to  7  digits  at  the  termination,  where  the  maximum  allowed  iterations  per  step  in 
DASSL  has  been  increased  to  50.  In  Table  2,  using  initial  stepsize  ho  =  5.0  x  10“3, 
starting  values  x0  =  0.8,  y0  =  -yft  -  *o  =  0,  and  zero  velocity  vector,  the  CS 
and  CM  iterations  are  nearly  identical  as  the  solution  TOL  decreases  from  10-4  to 
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MethocT 

ho 

no.  iter. 

approx,  rate  oj  convg. 

ur~ 

10-J 

24 

0.31602489817881 

CM 

10“3 

24 

0.31602542368461 

cs 

5  x  10~3 

46 

0.52240576939833 

CM 

5  x  10"3 

46 

0.52240703666124 

CS 

io-2 

74 

0.65889119756221 

CM 

io-2 

74 

0.65889222379467 

CS 

5  x  IO"2 

96 

0.72329054568167 

CM 

5  x  IO"2 

96 

0.72329103128307 

Table  1:  Simplified  2D  Bushing  Force  -  Rates  of  Convergence  for  Nonlinear  Iteration 
with  Initial  Stepsize  ho 


MethocT 

TOE" 

step 

j~s 

etj  —  s 

7tT=T 

CS 

io-4 

17 

79 

17 

3 

1 

CM 

IO"4 

17 

79 

17 

3 

1 

CS 

10“6 

17 

130 

17 

2 

2 

CM 

10~6 

17 

129 

17 

2 

2 

CS 

IO"8 

13 

146 

11 

0 

2 

CM 

IO"8 

13 

147 

11 

0 

2 

CS 

10-i° 

13 

175 

11 

0 

2 

CM 

10-1° 

13 

176 

11 

0 

2 

CS 

10-12 

13 

200 

11 

0 

2 

CM 

10-12 

13 

202 

11 

0 

2 

Table  2:  Simplified  2D  Bushing  Force  -  Results  from  DASSL  in  0-0.05  Seconds  Run 
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10-12  in  the  simplified  bushing  example  using  the  first-order  Euler  method  in  DASSL. 
Because  of  the  strong  damping  property  of  the  backward  Euler  method,  the  solution 
converged  to  the  fixed  point  (1,  — 10-5, 0)  after  the  time  t  =  0.05  seconds. 

We  also  examined  the  case  where  the  bushing  force  continuously  generates  oscilla¬ 
tory  forces  so  that  the  predictor  may  be  far  away  from  the  smooth  solution.  Preloading 
the  bushing  force  by  moving  the  attachment  points  from  (|,0)  and  (— §,0)  to  (0,0) 
on  the  global  and  the  body-fixed  frames,  respectively,  i.e.,  the  bushing  force  becomes 


r/xi 

’  x  ' 

— m 

y 

j9 . 

$ 

For  q  =  (x,  y,  6)  satisfying  the  constraints,  the  magnitude  of  the  bushing  force  is  O(^), 
since  6=  ||?— ?*||  =  0(1).  Using  the  initial  values  [1,0, 0,0, 0,0],  /*  =  1  ,TOL  =  10“10, 
and  e  =  10“5,  Tab.  3  contains  the  results  of  the  CS  and  CM  iterations.  To  observe 
the  behavior  of  Newton  convergence,  we  have  chosen  the  stepsize  ho  from  10-3  to 
5  x  10-2  and  ending  time  equal  to  h0,  with  the  order  of  BDF  in  DASSL  restricted  to 
one.  For  this  test,  DASSL  was  modified  so  that  the  number  of  Newton  iteration  per 
time  step  is  <  100  (instead  of  the  default  4). 

Under  these  conditions,  the  CS  and  CM  iterations  are  recorded  in  Tab.  3.  The 
rate  of  convergence  of  the  CM  iteration  is  commensurable  to  the  stepsize,  i.e.,  the 
number  of  function  evaluations  decreases  when  the  stepsize  is  halved.  For  larger  step- 
sizes,  i.e.,  ho  ranging  from  0.005  to  0.00075,  the  CS  iteration  exhibits  more  frequent 
Jacobian  changes  than  the  CM  iteration  because  of  the  error  and  convergence  test 
failures  in  the  CS  iterations.  Hence,  more  than  one  steps  were  taken  when  DASSL 
restarted  with  a  smaller  stepsize  This  agrees  with  the  prediction  given  by  Theo¬ 
rem  3,  since  the  constrained  potential  force  is  of  the  same  order  of  magnitude  as  the 
bushing  force,  i.e.,  0(|).  The  CS  iteration  performs  much  better  for  smaller  hQ,  e.g., 
from  0.0005  to  10"4,  than  for  larger  h0.  This  is  because  for  h  <  >/e,  \0(h2)  «  0(h2) 
for  the  forward  Euler  predictor.  For  the  case  ho  =  0.0001,  the  CS  iteration  converged 
in  one  iteration. 

In  a  full  2D  bushing  example,  i.e.,  the  two-body  bi-directional  force  (44),  which 
was  preloaded  as  the  simplified  bushing,  a  large  number  of  convergence  test  failures 
in  the  CS  iterations  was  observed  in  Table  4.  For  this  two-body  model,  we  used 
g  =  13.5,  ll  =  1,  **  =  kv  =  ke  =  10-5,  and  c*  =  c*  =  c9  =  10.  Using  the 
initial  values  [0, 0, 0, 0, 0.9989,  —0.001485, 0]  and  the  initial  velocity  [0, 0, 0,  — 0.675e  — 
4,  — 0.45444e  — 2,0],  results  from  DASSL  of  a  0  to  0.5  second  simulation  are  contained 
in  Table  4.  In  this  case,  even  when  the  position  and  velocity  are  within  the  TOL  as 
shown  in  the  last  two  columns  in  Tab.  3,  |  can  be  large  since  the  preloaded  bushing 
may  cause  ^  >  t. 
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I  Method" 

h~ 

E&3 

im 

J~s 

ofi a 

Kfl 

W£fl 

2 

i 

CM 

0.005 

H 

57 

1 

0 

0 

CS 

0.0025 

u 

80 

9 

2 

i 

CM 

0.0025 

u 

20 

1 

0 

0 

CS 

0.001 

mB 

28 

7 

2 

0 

CM 

0.001 

i 

8 

1 

0 

0 

CS 

0.00075 

9 

52 

12 

5 

0 

CM 

0.00075 

1 

6 

1 

0 

0 

CS 

0.0005 

8 

26 

4 

2 

0 

CM 

0.0005 

1 

4 

1 

1  0 

0 

Table  3:  Preloaded  Simplified  2D  Bushing  Force 


Method 

tut 

step 

7^ 

J  -  * 

etf  -  s 

ctf  -s 

X 

V 

CS 

~W3~ 

125 

465 

147 

19 

28 

0.8531799 

-0.6712374 

CM 

io-3 

27 

55 

9 

0 

0 

0.8527 

-0.6431779 

CS 

io-“ 

96 

278 

86 

2 

22 

0.8528879 

-0.6433484 

CM 

io-4 

35 

72 

11 

1 

0 

0.8527338 

-0.6447480 

CS 

10"5 

188 

630 

225 

1 

57 

0.8527399 

-0.6447200 

CM 

10"5 

45 

94 

7 

0 

0 

0.8527340 

-0.6447466 

CS 

10"6 

172 

542 

148 

2 

39 

0.8527345 

-0.6447471 

CM 

io-6 

60 

122 

8 

1 

0 

0.8527342 

-0.6447488 

Table  4:  Two-body  2D  Bushing  Force  in  0-0.5  Seconds  Run 
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